Libraries & Functions

Let’s begin by loading the libraries and functions that will be used for this analysis:

##                      plyr                     dplyr 
##                      TRUE                      TRUE 
##                data.table                    readxl 
##                      TRUE                      TRUE 
##                   stringr                   stringi 
##                      TRUE                      TRUE 
##                  forecast                 tidyverse 
##                      TRUE                      TRUE 
##                     caret                       zoo 
##                      TRUE                      TRUE 
##                    plotly                        DT 
##                      TRUE                      TRUE 
##                     rpart                     tidyr 
##                      TRUE                      TRUE 
##                 inspectdf              DataExplorer 
##                      TRUE                      TRUE 
##                 ggfortify                  magrittr 
##                      TRUE                      TRUE 
##                   ggpmisc                   prophet 
##                      TRUE                      TRUE 
##                   ggplot2                tsoutliers 
##                      TRUE                      TRUE 
##                       dlm                  reshape2 
##                      TRUE                      TRUE 
## AppliedPredictiveModeling      PerformanceAnalytics 
##                      TRUE                      TRUE 
##                  ggdendro                     rgdal 
##                      TRUE                      TRUE 
##                   tsibble 
##                      TRUE

Data Collection

Raw data for each Peru COVID-19 dataset are loaded.

#load raw data
## Load data
load("rda/COVID_19_Peru_Raw.rda")

load("rda/COVID_19_Peru_Raw_Inc2day.rda")

Data Summary

The COVID-19 dataset (COVID_19_Peru_Raw) has 39 rows and 12 columns while COVID-19 dataset case increase every 2 days has 20 rows and 3 columns

Exploratory Data Analysis (EDA)

After raw data is loaded, I proceed to perform EDA on dataset.

introduce(COVID_19_Peru_Raw)

## # A tibble: 1 x 9
##    rows columns discrete_columns continuous_colu~ all_missing_col~
##   <int>   <int>            <int>            <int>            <int>
## 1    39      12                1               11                0
## # ... with 4 more variables: total_missing_values <int>,
## #   complete_rows <int>, total_observations <int>, memory_usage <dbl>

introduce(COVID_19_Peru_Raw_Inc2day)

## # A tibble: 1 x 9
##    rows columns discrete_columns continuous_colu~ all_missing_col~
##   <int>   <int>            <int>            <int>            <int>
## 1    20       3                1                2                0
## # ... with 4 more variables: total_missing_values <int>,
## #   complete_rows <int>, total_observations <int>, memory_usage <dbl>

COVID_19_Peru_Raw$Date <- as.Date(COVID_19_Peru_Raw$Date, format = "dd/MM/YY")
COVID_19_Peru_Raw_Inc2day$Date <- as.Date(COVID_19_Peru_Raw_Inc2day$Date, format = "dd/MM/YY")

Peru_Daily_Incident <- COVID_19_Peru_Raw %>% pivot_longer(-Date, names_to = "province", 
    values_to = "incident_cases") %>% filter(province %in% c("Lima_cases_increase", 
    "Other_province_cases_increase", "Peru_cases_increase")) %>% mutate(province = ordered(province, 
    c("Lima_cases_increase", "Other_province_cases_increase", "Peru_cases_increase")))

Peru_Cumulative_Incidence <- Peru_Daily_Incident %>% group_by(province) %>% 
    arrange(Date) %>% tidyr::replace_na(list(incident_cases = 0)) %>% mutate(cumulative_incident_cases = cumsum(incident_cases)) %>% 
    mutate(cumulative_incident_cases = ifelse(cumulative_incident_cases == 0, 
        NA, cumulative_incident_cases))

# Plot graph
a <- Peru_Cumulative_Incidence %>% ggplot(aes(x = Date, y = cumulative_incident_cases)) + 
    geom_point(color = "darkblue") + geom_line() + scale_x_date(date_breaks = "7 days", 
    date_labels = "%d %b") + facet_grid(province ~ ., scales = "free_y") + labs(y = "Daily cumulative incidence", 
    title = "Lab-confirmed possitive cases of COVID-19 in Peru, 2020", caption = "Note: varying y-axis scales") + 
    theme(legend.position = "none", strip.text.y = element_text(size = 11)) + 
    theme_bw()

ggplotly(a)

b <- Peru_Cumulative_Incidence %>% ggplot(aes(x = Date, y = cumulative_incident_cases)) + 
    geom_point(color = "darkblue") + geom_line() + scale_y_log10() + scale_x_date(date_breaks = "7 days", 
    date_labels = "%d %b") + facet_grid(province ~ ., scales = "free_y") + labs(y = "Daily cumulative incidence (log10)", 
    title = "Lab-confirmed possitive cases of COVID-19 in Peru, 2020", caption = "Note: varying y-axis scales") + 
    theme(legend.position = "none", strip.text.y = element_text(size = 11)) + 
    theme_bw()

ggplotly(b)

c <- Peru_Daily_Incident %>%
  ggplot(aes(x=Date, y=incident_cases))  + # geom_point() + geom_line() +
    geom_bar(stat="identity", fill = "darkblue") + 
    scale_x_date(date_breaks="7 days", date_labels = "%d %b") +
    facet_grid(province ~., scales="free_y") + labs(y="Daily incremental incidence",
                                   title="Lab-confirmed possitive cases of COVID-19 in Peru",
                                   caption="Note: varying y-axis scales")  +
    theme(legend.position = "none", 
          strip.text.y = element_text(size=11)) + theme_bw()

ggplotly(c)

x <- COVID_19_Peru_Raw %>%
  ggplot(aes(x=Date, y=increase_index_case_percent_Peru))  + # geom_point() + geom_line() +
    geom_bar(stat="identity", fill = ifelse(COVID_19_Peru_Raw$increase_index_case_percent_Peru < 0, "darkred",
                                            "darkblue")) + 
    scale_x_date(date_breaks="7 days", date_labels = "%d %b") +
    labs(y="% Index of possitive cases per day",
         title="Percentage Index for Lab-confirmed possitive cases of COVID-19 in Peru")  +
    theme(legend.position = "none", 
          strip.text.y = element_text(size=11)) + theme_bw() +
  scale_y_continuous(labels = scales::percent)

ggplotly(x)

z <- COVID_19_Peru_Raw_Inc2day %>%
  ggplot(aes(x=Date, y=percent_cases_increase_2_days))  + 
  geom_point(color = "darkblue") + geom_line() +
    #geom_bar(stat="identity", fill = "darkblue") + 
    scale_x_date(date_breaks="2 days", date_labels = "%d %b") +
    labs(y="Incremental incidence every 2 days",
         title="Percentage of Lab-confirmed possitive cases of COVID-19 in Peru every 2 days")  +
    theme(legend.position = "none", 
          strip.text.y = element_text(size=11)) + theme_bw() +
  scale_y_continuous(labels = scales::percent)

ggplotly(z)

# Change the 0 on deaths to NA

COVID_19_Peru_Raw <- COVID_19_Peru_Raw %>% mutate(Deaths_Total = ifelse(Deaths_Total == 
    0, NA, Deaths_Total)) %>% mutate(Death_Daily = ifelse(Death_Daily == 0 & 
    Date < ymd("2020-03-19"), NA, Death_Daily)) %>% mutate(percent_death_increase_per_day = ifelse(percent_death_increase_per_day == 
    0 & Date < ymd("2020-03-19"), NA, percent_death_increase_per_day)) %>% mutate(`increase_index_death_percent)Peru` = ifelse(`increase_index_death_percent)Peru` == 
    0 & Date < ymd("2020-03-19"), NA, `increase_index_death_percent)Peru`))


d <- COVID_19_Peru_Raw %>% ggplot(aes(x = Date, y = Deaths_Total)) + geom_bar(stat = "identity", 
    fill = "brown") + scale_x_date(date_labels = "%d %b") + labs(y = "Daily cumulative deaths", 
    title = "Cumulative deaths in lab-confirmed possitive cases of COVID-19 in Peru", 
    caption = "Note: varying y-axis scales") + theme(legend.position = "none", 
    strip.text.y = element_text(size = 11)) + theme_bw()

ggplotly(d)

e <- COVID_19_Peru_Raw %>% ggplot(aes(x = Date, y = Death_Daily)) + geom_bar(stat = "identity", 
    fill = "brown") + scale_x_date(date_labels = "%d %b") + labs(y = "Daily incremental deaths", 
    title = "Daily deaths in lab-confirmed possitive cases of COVID-19 in Peru", 
    caption = "Note: varying y-axis scales") + theme(legend.position = "none", 
    strip.text.y = element_text(size = 11)) + theme_bw()

ggplotly(e)

y <- COVID_19_Peru_Raw %>% mutate(colorbars = ifelse(COVID_19_Peru_Raw$`increase_index_death_percent)Peru` > 
    0, "darkred", "darkblue")) %>% ggplot(aes(x = Date, y = `increase_index_death_percent)Peru`, 
    fill = colorbars)) + geom_bar(stat = "identity") + theme_bw() + scale_x_date(date_breaks = "7 days", 
    date_labels = "%d %b") + labs(y = "% Index of deaths per day", title = "Percentage Index for Death Lab-confirmed possitive cases of COVID-19 in Peru") + 
    theme(legend.position = "none") + scale_y_continuous(labels = scales::percent)

ggplotly(y)

Using the SIR model for Peru to make predictions

An obvious next step is to use our fitted SIR model to make predictions about the future course of the outbreak. However, caution is required, because the SIR model assumes a fixed reproduction number, but if public health interventions have been implemented, such as quarantining of cases, contact tracing and isolation of those contacts, and general restrictions on social mixing, such as closing the country, then the effective reproduction number \(R_{e}\) will be dynamic and should fall as those interventions are progressively implemented, to values considerably less than the basic reproduction number \(R_{0}\), which reflects the behaviour of the virus at the beginning of an epidemic before any response has been implemented.

So let’s use our SIR model, fitted to the first 39 days of data, to extrapolate out to the current date, and compare that against the observed values:

We can see that the actual incidence similar than that predicted by our model. The reason is that, even with public health interventions implemented by the Peruvian authorities, the \(R_{e}\) of the COVID-19 in Peru hasn’t been reduced by much. This is something that will need to be addressed ASAP.

When the \(R_{e}\) falls below 1.0, the peak of the epidemic will have been reached and the outbreak will eventually die out.

Using our model to let the outbreak “run its course” without intervention

It is instructive to use our model fitted to the first 38 days of available data on lab-confirmed cases in Peru, to see what would happen if the outbreak were left to run its course, without public health interventions.

# time in days for predictions
t <- 1:100
# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(y = init, times = t, func = SIR, 
    parms = Opt_par))
# add a Date column and join the observed incidence data
fitted_cumulative_incidence <- fitted_cumulative_incidence %>% mutate(Date = ymd(sir_start_date) + 
    days(t - 1), province = "Peru_cases_increase") %>% left_join(Peru_Cumulative_Incidence %>% 
    ungroup() %>% filter(province == "Peru_cases_increase") %>% select(Date, 
    cumulative_incident_cases))
# plot the data
h <- fitted_cumulative_incidence %>% ggplot(aes(x = Date)) + geom_line(aes(y = I), 
    colour = "red") + geom_line(aes(y = S), colour = "black") + geom_line(aes(y = R), 
    colour = "green") + geom_point(aes(y = cumulative_incident_cases), colour = "blue") + 
    scale_y_continuous(labels = scales::comma) + labs(y = "Persons", title = "COVID-19 fitted vs observed cumulative incidence, Peru", 
    subtitle = "(red = Infectious, black = Susceptible, green = Recovered, blue = Observed incidence") + 
    scale_colour_manual(name = "", values = c(red = "red", black = "black", 
        green = "green", blue = "blue"), labels = c("Susceptible", "Recovered", 
        "Observed incidence", "Infectious")) + theme_bw()

ggplotly(h)

It is easier to see what is going on if we use a log scale:

Clearly that prediction, should it come to pass, would be an unmitigated disaster. At this point, it is worth remarking on the importance of decisive public health intervention to limit the spread of such epidemics. Without such interventions, tens of millions of people could be infected, as our model predicts, and even with only a one or two per cent mortality rate, hundreds of thousands of people would die. The economic and human cost of the outbreak control steps taken within China are large, but the alternatives are a lot worse!

Ascertainment rates

So far, we have assumed that the counts of lab-confirmed cases represent all the cases that are infectious. This is unlikely to be true – typically only a proportion of actual cases are detected or found or sent for testing. This proportion is known as the ascertainment rate. The ascertainment rate is likely to change during the course of an outbreak, particularly if surveillance and screening efforts are increased, or if case definitions are changed. Such changing ascertainment rates can be easily incorporated into our model by using a set of weights, or a weighting function, for our incidence data, but for the sake of simplicity, let’s see what happens if we assume a fixed ascertainment rate of 20%². If we apply that, thus inflating the number of incident cases by a factor of 5, and refit our model, we get the following results.

# put the daily cumulative incidence numbers for Hubei from 6th Mar to 14th
# Apr into a vector called Infected
sir_start_date <- "2020-03-07"
Infected <- Peru_Cumulative_Incidence %>% filter(province == "Peru_cases_increase", 
    Date >= ymd("2020-03-07"), Date <= ymd("2020-04-14")) %>% pull(cumulative_incident_cases)
# Apply a fixed 20% ascertainment rate
Infected <- Infected * 5
# Create an incrementing Day vector the same length as our cases vector
Day <- 1:(length(Infected))
# now specify initial values for S, I and R
init <- c(S = N - Infected[1], I = Infected[1], R = 0)
RSS <- function(parameters) {
    names(parameters) <- c("beta", "gamma")
    out <- ode(y = init, times = Day, func = SIR, parms = parameters)
    fit <- out[, 3]
    sum((Infected - fit)^2)
}
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 
    1))
# check for convergence
Opt$message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par

##      beta     gamma 
## 0.6031149 0.3968851

R0 <- setNames(Opt_par["beta"]/Opt_par["gamma"], "R0")
R0

##       R0 
## 1.519621

Note that these fitted parameters are the same as the ones we got above, without an ascertainment rate adjustment. Let’s look at the fitted values.

# time in days for predictions
t <- 1:100
# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(y = init, times = t, func = SIR, 
    parms = Opt_par))
# add a Date column and join the observed incidence data
fitted_cumulative_incidence <- fitted_cumulative_incidence %>% mutate(Date = ymd(sir_start_date) + 
    days(t - 1), province = "Peru_cases_increase") %>% left_join(Peru_Cumulative_Incidence %>% 
    ungroup() %>% filter(province == "Peru_cases_increase") %>% select(Date, 
    cumulative_incident_cases))
# plot the data

i <- fitted_cumulative_incidence %>% ggplot(aes(x = Date)) + geom_line(aes(y = I), 
    colour = "red") + geom_line(aes(y = S), colour = "black") + geom_line(aes(y = R), 
    colour = "green") + geom_point(aes(y = cumulative_incident_cases * 5), colour = "blue") + 
    scale_y_log10(labels = scales::comma) + labs(y = "Persons", title = "COVID-19 fitted vs observed cumulative incidence adjusted for 20% ascertainment, Peru") + 
    scale_colour_manual(name = "", values = c(red = "red", black = "black", 
        green = "green", blue = "blue"), labels = c("Susceptible", "Recovered", 
        "Observed incidence", "Infectious")) + theme_bw()

ggplotly(i)

Perhaps counter-intuitively, incorporation of a fixed 20% ascertainment rate adjustment into our model makes little difference to the modelled outbreak if it is let run its course, except that it all happens a bit more quickly. But the number of infections remains the same, and the basic reproduction number is unchanged. Note that that is for a fixed ascertainment rate. If the ascertainment rate varies significantly over time, then the parameter estimates will necessarily be biased – but in the early days of an outbreak, it may be reasonable to assume that ascertainment rates don’t change too much.